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ABSTRACT 


We  are  developing  a  new  algorithm  for  seismic  event  location  to  address  particular  problems  of 
importance  in  monitoring  the  Comprehensive  Nuclear-Test-Ban  Treaty,  including  the  assessment 
of  location  accuracy  and  increasing  the  number  of  seismic  phases  used  for  location.  Our  approach 
is  formulated  within  a  maximum-likelihood  estimation  framework  and  implemented  numerically 
with  grid-search  and  Monte  Carlo  techniques.  It  obtains  globally  optimal  hypocentral  estimates 
and  non-ellipsoidal  confidence  regions  that  do  not  depend  on  the  usual  assumption  of  local  linearity 
of  the  forward  problem.  The  approach  accommodates  such  complexities  as  3-D  traveltime  tables, 
non-Gaussian  data  errors,  and  general  types  of  parameter  constraints.  Specific  problems  we  are 
addressing  are  (1)  extension  of  our  uncertainty  analysis  to  account  more  realistically  for  the  effects 
of  modeling  errors  (errors  in  traveltime  tables),  which  traditionally  have  not  been  distinguished 
from  random  picking  errors  in  the  determination  of  location  confidence  regions;  (2)  incorporation 
of  azimuth  and  slowness  data,  in  addition  to  arrival  times,  into  the  maximum-likelihood  framework; 
and  (3)  integration  of  phase  association  into  the  event  location  algorithm  in  order  to  perform  these 
tasks  simultaneously  with  a  grid-search  approach.  We  present  examples  of  our  event  location 
algorithm  applied  to  data  from  the  International  Monitoring  System,  illustrating  the  effects  of 
nonlinearity  on  hypocentral  confidence  regions  determined  from  sparse  data  sets. 
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OBJECTIVE 


Over  the  past  decade,  grid  search  has  gained  popularity  as  an  inversion  algorithm  for  seis¬ 
mic  event  location,  owing  to  its  simplicity  of  implementation  and  ability  to  find  globally  optimal 
hypocentral  estimates.  Unlike  the  traditional  iterative  methods,  grid  search  easily  accommodates 
problem  complexities  such  as  general  error  models  (non-Gaussian),  nonlinear  parameter  constraints, 
and  multiple  data  types  (azimuth  and  slowness  in  addition  to  arrival  times).  One  of  the  early  ap¬ 
plications  of  the  method  was  the  work  of  Sambridge  and  Kennett  (1986),  who  applied  a  directed 
grid-search  technique  to  teleseismic  location.  Recently,  Dreger  et  al.  (1998)  applied  grid  search  to 
problems  of  event  location  with  sparse  regional  networks,  an  important  situation  in  CTBT  moni¬ 
toring.  Previous  applications  of  grid  search  in  event  location  were  formulated  to  minimize  a  “data 
misfit”  function  with  respect  to  hypocentral  parameters.  Contours  of  this  function  were  recognized 
to  approximate  the  boundaries  of  confidence  regions,  to  which  confidence  levels  could  be  assigned 
within  certain  linear  approximations  and  under  the  assumption  of  Gaussian  errors  (Wilcock  and 
Toomey,  1991). 
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We  are  developing  a  more  general  algorithm  for  event  location  that  uses  grid  search  to  maximize 
a  general  likelihood  function.  Maximum  likelihood  estimation  theory  provides  a  formal  framework 
for  non-Gaussian  error  assumptions  and  for  combining  diverse  types  of  data  and  prior  information 
on  parameters.  In  addition,  our  algorithm  replaces  approximate,  analytical  formulas  for  confidence 
levels  with  Monte  Carlo  simulation,  which  is  not  restricted  by  the  Gaussian  assumption  and  which 
accounts  fully  for  nonlinearity  of  the  forward  problem  and  parameter  constraints.  Other  enhance¬ 
ments  in  progress  or  planned  for  the  near  future  include  the  incorporation  of  azimuth  and  slowness 
data,  a  more  complete  treatment  of  the  errors  in  traveltime  tables  sued  for  forward  modeling,  and 
integration  of  seismic  phase  association  tasks  into  our  grid-search  algorithm  in  an  attempt  to  con¬ 
strain  event  locations  with  more  seismic  phases.  Our  objective  is  improved  location  accuracy  and 
more  realistic  analyses  of  location  error  for  small,  sparsely  recorded  events. 

RESEARCH  ACCOMPLISHED 


Maximum  Likelihood  Formulation 


The  hypocentral  parameters  of  a  seismic  event  are  a  three-dimensional  position  vector  x  and  an 
origin  time  t.  Let  d  =  (di,  cfo-  ■  ■  ■ ,  dn)  be  an  n-dimensional  vector  of  arrival  times  picked  from 
various  seismic  phases  at  a  seismic  network.  The  event  location  problem  may  be  expressed  as 


<k  =  t  +  Tj(x)  +  e*,  *  =  l,...,n 


(1) 


where  T)  is  a  traveltime  function  (traveltime  table)  for  the  ith  datum  and  e,  is  an  error  with  respect 
to  this  function.  The  index  i  counts  over  both  stations  and  phase  types  (P,  S,  etc.),  including  only 
those  combinations  that  have  been  actually  observed  (i.e. ,  completeness  of  data  coverage  is  not 
assumed) . 

We  assume  the  errors  are  statistically  independent  and,  following  Billings  et  al.  (1994),  that 
each  is  distributed  with  a  “generalized  Gaussian”  probability  density  function  (p.d.f.)  given  by 
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where  K  =  2p1/JT(l  +  1/p)  and  Y  is  the  gamma  function.  When  p  =  2,  then  e%  is  normally 
distributed  with  zero  mean  and  variance  (crj)2.  When  p  —  1,  it  is  exponentially  distributed.  We 
assume  that  the  scale  parameters  at  are  known  in  a  relative  sense  and  write 


Ui  =  avi 


(3) 


where  the  //,;  are  known  but  the  universal  scale  parameter,  <r,  is  not. 

The  joint  p.d.f.  of  the  n  data  is  the  product  of  the  error  p.d.f.’s.  Considered  as  a  function  of 
the  unknown  parameters  (x,  t  and  a)  this  joint  p.d.f.  is  a  likelihood  function,  which  we  denote  as 
L(x,i,  <r;d).  It  is  convenient  to  deal  with  the  negative  logarithm  of  likelihood,  which  we  denote  as 
A.  This  is  given  by 


A(x,  t.  cr;  d) 


-  log  L(x,  t,  <j;  d) 
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Here,  4/  is  a  ‘data  misfit’  function  defined  as 
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The  maximum  likelihood  estimates  of  the  unknowns  are  the  values  that  maximize  L  or,  equivalently, 
minimize  A.  We  denote  these  estimates  as  xmj,  tmi,  and  am[.  The  maximization  may  be  subjected 
to  prior  constraints  on  the  parameters.  In  our  current  algorithm,  we  assume  hard  upper  and  lower 
bounds  on  focal  depth  {z)  and  the  p.d.f.  scale  factor: 

0  <  z  <  zmax 
amin  <o<  amax 

We  allow  latitude  and  longitude  to  be  anywhere  on  the  globe,  or  restricted  to  within  a  given  distance 
of  a  given  geographic  location.  We  place  no  constraint  on  origin  time. 

Given  its  structure,  A  is  amenable  to  a  hierarchical  minimization  with  respect  to  the  unknown 
parameters,  leading  to  minimization  subproblems  of  lower  dimension  than  the  original  problem.  We 
define  a  “reduced”  objective  function  which,  for  each  fixed  hypocenter,  is  minimum  with  respect 
to  t  and  o  (subject  to  prior  constraints  on  a): 

A(x;  d)  =  min  A(x,  i,  er;  d).  (8) 

i,< 7 

The  location  problem  reduces  to  minimization  of  A. 


(6) 

(7) 


Grid  Search 

Our  current  grid  search  algorithm  obtains  the  maximum  likelihood  estimates  of  the  hypocentral 
parameters  (xmj  and  tmi)  and  the  p.d.f.  scale  parameter  ( ami ),  as  defined  in  the  previous  section. 
The  algorithm  computes  the  reduced  objective  function,  A  in  equation  (8),  at  each  point  in  a  3-D 
grid  of  hypocenters.  Following  previous  workers,  the  hypocenter  grid  is  constructed  dynamically 
through  a  process  of  successive  refinement.  Our  procedure  for  grid  refinement  resembles  that  of  the 
‘neighborhood’  search  algorithm  developed  by  Sambridge  (1999).  The  first  grid  covers  the  entire 
globe  and  from  0  to  700  km  in  depth  at  coarse  spacing:  100  km  in  depth,  9  degrees  in  latitude,  and 
9  degrees  in  longitude  near  the  equator  and  increasing  at  higher  latitudes.  On  each  pass  of  grid 
refinement,  nodes  are  added  as  neighbors  of  a  subset  of  grid  points  comprising  the  “best”  (smallest 
A)  points  tested  thus  far.  Neighbors  are  placed  at  one-third  the  grid-spacing  of  the  previous  pass. 
The  size  of  the  grid  subset  chosen  for  refinement  is  reduced  on  each  pass.  The  search  ends  when 
the  grid  spacing  is  less  than  0.3  km. 


Non-Ellipsoidal  Confidence  Regions 

A  confidence  region  on  the  hypocenter  of  an  event  is  defined  classically  in  terms  of  a  statistic,  r, 
that  is  a  function  of  the  data  vector,  d,  and  the  hypocenter,  x.  A  confidence  region  at  confidence 
level  1  —  a  comprises  those  values  of  x  satisfying  the  inequality 

cdf(r(d,  x))  <  1  —  a.  (9) 


Here,  cdf  denotes  the  cumulative  distribution  function  of  a  random  variable. 

Flinn  (1965)  defined  hypocentral  confidence  regions  for  the  case  of  Gaussian  errors  ( p  =  2)  and 
no  prior  bounds  on  a  (amm  =  0,  anrdx  =  oo).  Evernden  (1969)  treated  the  case  with  a  completely 
known  (a  =  amm  =  cjmdx  =  omi),  while  Jordan  and  Sverdrup  (1981)  generalized  these  two  results 
to  account  for  partial  prior  knowledge  of  a.  The  statistic  used  by  each  of  these  workers  can  be 
expressed  in  our  notation  as 

T(x;d)  -  T(xmJ;d) 


where  T  is  the  “reduced”  data  misfit  function  (T  minimized  with  respect  to  t  for  each  fixed  x). 
Under  the  approximation  that  the  traveltime  functions,  T*(x),  are  locally  linear  near  x  =  xmj,  the 
probability  distribution  of  r  does  not  depend  on  any  of  the  unknown  parameters  (x,  t  or  cr).  and 
a  confidence  region  takes  the  form 

§(x;d)  -  fk(xmJ;d)  <  a2mlCa  (11) 


where  Ca  is  a  constant.  Under  the  linearity  assumption,  this  equation  defines  a  hyper-ellipsoid  in 
hypocenter  space.  Flinn  (1965),  Evernden  (1969),  and  Jordan  and  Sverdrup  (1981)  derive  different 
values  of  Ca  based  on  differing  assumptions  about  a. 

To  test  the  linear  approximation,  Wilcock  and  Toomey  (1991)  computed  non-ellipsoidal  confi¬ 
dence  regions  directly  from  equation  (11).  This  equation  implies  that  each  surface  of  constant  data 
misfit,  T(x;d),  is  the  boundary  of  a  confidence  region  for  some  level  of  confidence,  1  —  a.  They 
directly  sampled  T(x;  d)  on  a  grid  of  hypocenters  to  map  the  non-ellipsoidal  confidence  region  for 
a  fixed  confidence  level,  thus  taking  the  major  effect  of  nonlinearity  into  account.  However,  the 
approach  of  Wilcock  and  Toomey  (1991)  ignores  another  effect  of  nonlinearity:  the  probability 
distribution  of  r  may  depend  on  the  true  values  of  the  unknown  parameters.  This  effect  can  be 
caused  by  the  use  of  nonlinear  constraints  on  the  parameters,  as  in  eqs.  (6)  and  (7),  as  well  as  by 
nonlinearity  of  the  traveltime  functions.  Moreover,  they  did  not  consider  the  case  of  non-Gaussian 
errors. 

We  have  generalized  the  approach  of  Wilcock  and  Toomey  as  follows.  We  define  a  test  statistic 
as  the  logarithm  of  a  likelihood  ratio: 


r(d,  x) 


log 


maxX;t,g  T(x,  t,  cr;  d) 
max^o-  L(x,  t,  cr,  d) 


A(x;  d)  -  A(xmJ;d). 


(12) 


That  is,  r  differences  our  reduced  objective  function  between  x  and  the  maximum  likelihood  solu¬ 
tion,  xmj.  This  statistic  is  equivalent  to  the  ones  used  by  Flinn  and  Evernden  under  their  respective 
assumptions  about  cr,  but  accommodates  arbitrary  constraints  on  the  parameters  and  non-Gaussian 
error  distributions.  Letting  x  =  (x,y,z),  the  statistic  for  a  2-D  confidence  region  on  the  event  epi¬ 
center,  (x,y),  is  given  by 

r(d,  x ,  y)  =  min  A(x;  d)  -  A(xmJ;  d)  (13) 

and  for  a  confidence  interval  on  focal  depth  is 

r(d,  z)  =  minA(x;d)  —  A(xmj;d).  (14) 

x>y 

Confidence  regions  using  these  log-likelihood  statistics  could  still  be  defined  via  the  inequality 
of  equation  (9),  except  this  inequality  presumes  that  the  distribution  (c.d.f.)  of  r  does  not  depend 
on  the  true  values  of  the  parameters.  We  assume  the  main  dependence  is  on  focal  depth  and  cr,  and 
write  the  c.d.f.  of  r  as  cdf(r;  z,  a).  We  generalize  the  inequality  of  equation  (9)  to  use  the  c.d.f.  of 
r  that  is  minimum  with  respect  to  the  true  parameters.  Thus,  the  hypocentral  confidence  region 
is  given  by 


mincdf(r(d, x);  £,  a)  <  1  —  a.  (15) 

(7 

The  epicentral  one  is 


mincdf(r(d, x, y)\ z, a)  <1  —a, 


(16) 


and  the  focal  depth  confidence  interval  is 


mincdf(r(d,  z)\  z,  a)  <  1  —  a.  (17) 

<r 

With  these  definitions,  a  confidence  region  will  include  the  true  value  of  the  parameter  at  least 
100(1  —  a)  percent  of  the  time. 

Confidence  Regions  Via  Monte  Carlo  Sampling 

We  outline  the  technique  for  hypocentral  confidence  regions.  The  basic  idea  is  to  estimate  the  c.d.f. 
of  the  test  statistic  r  by  simulation,  i.e.  computing  r  for  many  randomly  generated  samples  of  the 
error  vector.  We  generate  each  error,  ef*c,  using  a  pseudo-random  number  generator  in  accordance 
with  the  assumed  error  distribution  (eqs.  (2)-(3))  for  some  given  “true”  a.  Then,  for  given  true 
hypocentral  parameters,  x  and  t,  synthetic  data  are  calculated  as 

d™  =  t  +  Ti{x)  +  eflc.  (18) 

We  apply  our  grid  search  algorithm  to  these  data  to  obtain  the  m.l.  hypocenter,  x™j\  Plugging 
this  into  the  formula  for  r, 

T(dmc,x)  =  A(x;  dmc)  —  A(x™[;  dmc).  (19) 

we  obtain  one  sample  from  cdf (r;z,a).  We  compare  this  sample  to  the  observed  value  of  the 
statistic,  r(d,x),  obtained  from  the  real  data.  We  count  a  rejection  of  x  if 

r(<Tc,x)  <  r(d,x).  (20) 

The  fraction  of  rejections  after  many  Monte  Carlo  trials  yields  an  estimate  of  cdf(r(d,  x);  z,  a). 
Performing  this  simulation  for  multiple  values  of  a  and  then  minimizing  amongst  them  gives  the 
lowest  confidence  level  such  that  the  confidence  region  includes  x.  The  process  for  depth  confidence 
intervals  and  epicenter  confidence  regions  proceeds  in  the  same  manner,  except  that  in  the  latter 
case  the  simulation  is  performed  for  multiple  values  of  true  depth  as  well  as  a,  and  the  confidence 
level  is  minimized  over  both. 

Examples  With  IMS  Data 

We  illustrate  our  grid-search/Monte  Carlo  event  location  algorithm  on  two  events  from  the  Proto¬ 
type  International  Data  Center  (PIDC)  Reviewed  Event  Bulletin  (REB).  The  first  is  event  number 
20354875,  which  occurred  1  March  1999  near  Crete  and  was  assigned  a  magnitude  of  =  3.9. 
The  International  Monitoring  System  (IMS)  location  is  based  on  five  arrival  times:  Pn  at  stations 
BRAR  (A  =  7.2°)  and  KVAR  (A  =  15.7°),  and  P  at  BGCA  (A  =  30.6°),  ARCES  (A  =  34.5°) 
and  DBIC  (A  =  40.0°).  We  applied  our  grid-search  location  algorithm  to  this  event  using  the 
IASPEI91  traveltime  tables  and  with  no  elevation  or  path  corrections  applied  to  the  data.  Arrival 
time  errors  were  assumed  to  be  Gaussian  ( p  =  2)  with  standard  deviation  (cr)  bounded  between 
0.5  and  1.5  seconds. 

Figure  1  displays  confidence  level  as  a  function  of  event  location,  computed  via  Monte  Carlo 
simulation.  The  simulations  were  performed  with  three  hypothesized  “true”  values  of  a  and  6 
values  of  focal  depth.  The  number  of  error  realizations  for  each  a  and  z  was  300.  The  top  left  panel 
displays  confidence  level  as  a  function  of  event  depth.  The  two  points  where  the  curve  intersects  a 
given  confidence  level  are  the  end-points  of  the  focal  depth  confidence  interval  at  that  level.  The 


90%  confidence  interval  on  the  focal  depth  of  this  event  is  34  to  92  km,  while  the  95%  interval  is  0 
to  101  km.  The  IMS  solution  fixed  the  depth  to  zero. 

The  upper  right  panels  show  cross-sections  of  confidence  levels  on  the  3-D  hypocenter.  In  each 
plot,  a  contour  of  constant  confidence  level  outlines  the  boundary  of  the  hypocentral  confidence 
region  at  that  level.  The  REB  solution  for  the  hypocenter  is  plotted  as  a  circle.  We  see  that  the 
hypocentral  confidence  regions  display  significant  departure  from  ellipsoidal  behavior.  Significantly, 
the  effect  is  to  stretch  the  confidence  regions  towards  the  earth’s  surface,  making  it  more  difficult 
to  reject  a  shallow  focus. 

The  bottom  panels  are  confidence  level  as  a  function  of  event  epicenter,  whose  contours  define 
epicentral  confidence  regions  at  various  levels.  In  the  bottom  left  panel,  depth  is  unconstrained 
(z  >  0)  as  in  the  usual  definition  of  an  epicentral  confidence  region.  We  see  that  these  confidence 
regions  are  close  to  ellipsoidal  and  indicate  roughly  the  same  epicentral  uncertainty  as  the  REB, 
which  reports  semi-axes  of  57  and  22  km  for  the  90%  confidence  ellipse.  In  the  lower  right  panel, 
depth  is  constrained  to  be  shallow  (0  <  z  <  10  km).  The  title  above  the  plot  indicates  that  a 
shallow  focus  can  be  rejected  with  a  confidence  of  83%.  However,  the  plot  shows  that  this  low 
confidence  pertains  only  to  a  small  epicentral  region,  much  smaller  than  the  epicentral  confidence 
region  for  unconstrained  depth  (lower  left).  Thus,  the  epicenter  pertinent  to  a  shallow  focus  event 
is  much  better  determined. 

A  second  example  of  Monte  Carlo  confidence  regions  is  shown  in  Figure  2.  These  are  for  an 
rrif,  =  3.7  event  in  the  Afghanistan- Tajikistan  border  region  (17  March  1999,  event  no.  20380677 
from  PIDC).  The  data  set  comprised  five  P  wave  arrival  times,  all  teleseismic  except  station  NIL 
(A  =  5.1°).  The  REB  solution  reports  semi- axes  of  92  and  33  km  for  the  90%  confidence  region  on 
epicenter.  This  is  consistent  with  but  somewhat  larger  than  the  epicentral  uncertainty  implied  in 
the  lower  left  panel.  The  REB  confidence  interval  on  depth  is  35  to  130  km,  compared  to  our  smaller 
52  to  116  km  confidence  interval  (upper  left  panel  of  Figure  2).  Our  smaller  location  uncertainties 
may  be  due  to  a  more  optimistic  assumption  about  data  errors  (0.5  to  1.5  s)  than  used  in  the  IMS 
location  algorithm.  Looking  at  the  lower  right  panel,  we  see  that  a  shallow  focus  is  rejected  with 
98%  confidence,  and  essentially  with  100%  confidence  outside  an  elliptical  region  northeast  of  the 
unconstrained  solution. 

CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  developed  a  general  theoretical  and  computational  framework  for  finding  event  location 
estimates  and  confidence  regions,  relaxing  some  of  the  assumptions  currently  used  in  the  routine 
processing  of  IMS  data.  In  particular,  we  do  not  assume  local  linearity  of  forward  problem  (travel¬ 
time  vs.  hypocenter)  in  the  computation  of  confidence  regions,  and  we  accommodate  fairly  general 
error  models  and  parameter  constraints.  Our  preliminary  tests  show  consistency  with  IMS  results 
reported  in  the  Reviewed  Event  Bulletin,  but  also  some  noticeable  departures  from  the  ellipsoidal 
confidence  regions  computed  by  IMS  when  small  data  sets  are  used. 

Currently  we  are  incorporating  azimuth  and  slowness  data  into  our  formulation  and  algorithms. 
Such  data  are  routinely  used  in  IMS  processing  and  can  provide  important  constraints  on  sparsely 
recorded,  small  events.  We  are  beginning  to  address  some  additional,  difficult  problems  of  im¬ 
portance  in  CTBT  monitoring.  One  is  the  effect  of  errors  in  station  corrections  and  traveltime 
tables  (“modeling”  errors)  on  location  accuracy.  In  current  location  algorithms,  modeling  errors 
are  accounted  by  inflating  the  variances  assigned  to  the  data,  thus  assuming  they  are  random  and 
independent  between  stations  and  phases.  We  will  investigate  more  realistic  treatments  of  these 
errors  that  are  consistent  with  empirical  analyses  of  traveltime  residuals  and  calibration  studies. 
The  second  problem  to  address  is  unassociated  and  mis-associated  seismic  phases.  Our  goal  is  to 


integrate  key  elements  of  seismic  phase  association  into  the  process  of  event  location  in  an  attempt 
to  constrain  locations  with  more  data  and,  accordingly,  leave  fewer  unassociated  phases  in  the  ar¬ 
rival  bulletins.  Our  maximum  likelihood  formulations  and  grid-search  approach  to  event  location 
provides  a  flexible  framework  for  addressing  these  problems. 
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Figure  1:  Confidence  level  vs.  location  for  an  TO(,  =  3.9  event  near  Crete  (1  March  1999,  event 
no.  20354875  from  PIDC).  Top  left:  Confidence  level  vs.  focal  depth.  Top  center  and  right:  cross- 
sections  of  confidence  level  vs.  hypocenter.  Bottom  left:  Confidence  level  vs.  epicenter.  Bottom 
right:  Confidence  level  vs.  epicenter  with  focal  depth  constrained  between  0  and  10  km.  In  the  top 
images,  each  contour  of  constant  confidence  level  intersects  the  boundary  of  the  3-D  hypocentral 
confidence  region  at  that  level.  Contours  in  the  bottom  images  are  the  boundaries  of  epicentral 
confidence  regions.  Note  that  confidence  levels  below  80%  are  all  displayed  with  white.  The  circles 
mark  the  event  hypocenter  reported  in  the  Reviewed  Event  Bulletin. 
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Figure  2:  Confidence  level  vs.  location  for  an  =  3.7  event  in  the  Afghanistan- Tajikistan  border 
region  (17  March  1999,  event  no.  20380677  from  PIDC).  The  figure  is  in  the  same  format  as  Figure 
1.  Contours  of  constant  confidence  level  correspond  to  confidence  regions  on  the  hypocenter  (top 
images)  and  epicenter  (bottom  images).  The  circles  mark  the  event  hypocenter  reported  in  the 
Reviewed  Event  Bulletin. 


